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Abstract 

We developed a method to calculate positions and widths of three-body resonances. The method 
combines the hyperspherical adiabatic approach, slow variable discretization method (Tolstikhin 
et ai, J. Phys. B: At. Mol. Opt. Phys. 29, L389 (1996)), and a complex absorbing potential. 
The method can be used to obtain resonances having short-range or long-range wave functions. In 
particular, we have applied the method to obtain very shallow three-body Efimov resonances for a 
model system (Nielsen et ai, Phys. Rev. A 66, 012705 (2002)). 

PACS numbers: 31.15.Ja, 02.70.-c, 33.20.Tp 
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I. INTRODUCTION 



The problem of weakly bound states of three resonantly interacting particles, first dis- 
cussed long ago by Efimov has attracted a particular interest in recent years 
owing to the development of refined experimental techniques allowing tunning of the inter- 
particle interaction. Now, it has become possible to tune the two-body interaction in such a 
way that a two-body weakly bound or virtual state appears. In this situation, three particles 
interacting pairwise through the two-body potential have a series of weakly bound states 
If it happens that in addition to the weakly bound or virtual two-body state, there is 
another two-body bound state lying deeply below, the Efimov states are able to dissociate 
into a dimer and a free atom. Therefore, in such a situation, the Efimov states have finite 
lifetimes : these are Efimov resonances . Recently, such Efimov resonances have 

been observed in experiment 6]. 

n 

Efimov predicted an universal behavior for the energies of the Efimov states [1[ : If one 
knows the position, let's say, of the lowest state, then the others can be found using a scaling 
law 

£ n = £ n -iexp(27r/£), (1) 

where £ is a constant depending on the actual interaction potential. This is possible because 
the long-range potential of an Efimov-type system behaves as (£ 2 + 0.25)/(2/ip 2 ), where p 
is the hyper-radius and fi is the reduced three-body mass. The long-range behavior of the 
potential does not depend on details of the short-ranee two-body interaction. For the same 

un 

reason, widths of the Efimov resonances are scaled in the same way [7|, |8[ . Thus, if one needs 
to find positions and widths of the Efimov resonances for a particular three-body system, 
one has to be able to obtain at least one state: The approximate positions (and widths) of 
others can be obtained using the scaling law. 

Wave functions for such weakly bound Efimov states are characterized by a probability 
density located mainly at large inter-particle distances. On the other hand, in order to 
correctly obtain the accumulated phase, the small-distance part of the wave functions must 
also be represented accurately, even though the amplitude at such distances can be very 
small. Due to the very different lengthscales in the two regions, the calculation of the wave 
functions (of bound or resonant states) becomes difficult. 

Q 

Pen'kov 7] previously presented a theory of Efimov resonances using the formalism of 
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Faddeev equations. He calculated six resonances for a model two-body interaction potential: 
the obtained positions and widths are in agreement with the scaling law. Later, Nielsen et 
al. [8] computed Efimov resonances for a different model potential, using an adiabatic 
hyperspherical method coupled with a multi-channel i?-matrix method along the hyper- 
radius. In this approach, the authors obtained the scattering matrix as a function of energy 
above the lowest dissociation threshold. Then, they deduced widths and positions of the 
resonances from the variation of the phase shift 5(E) in the scattering matrix. 

Here, we are suggesting a different approach to calculate widths and positions of three- 
body resonances. This approach is general and can be used for weakly bound (long-range) 
and deeply-bound (short-range) resonances. The method uses hyperspherical coordinates 
similar to Ref. [8]. Instead of solving a multi-channel scattering problem, we suggest em- 
ploying the slow variable discretization (SVD) j^] method with adapted grid step [h]] in 
hyper-radius and hyperangles. As a basis set in hyper-radius we use the sine-DVR (dis- 
crete variable representation) method 11 , l^j]. At a large hyper-radius, we place a complex 



absorbing potential (CAP) to absorb the dissociating flux. Positions and widths of the 
three-body resonances are obtained as real and imaginary parts of complex eigenvalues of 
the resulting Hamiltonian of the system. 

The article is organized in the following way. In section [Til we present our approach of 
calculation of the three-body resonances. Then, in section IHIt we apply the method to 
obtain a single resonance of three particles interacting through a potential with a barrier 
[l^. The wave function of the resonance has a large amplitude mainly at relatively short 
distances. The second application of the developed method is discussed in section IIVI where 
we calculate several Efimov resonances for three particles interacting through a two-body 
short-range potential with a large two-body scattering length. We compare our results with 
a previous theoretical study [8|. Finally, in section [V] we discuss advantages and possible 
applications of the method. 
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II. METHOD 



A. Hyperspherical coordinates used in the present method 

In our approach we assume that the total angular momentum of the three-body system 
is zero. Thus, only three interparticle coordinates determine the configuration of the sys- 
tem. We represent the three interparticle coordinates by the Smith- Whitten hyperspherical 
coordinates [l4|. The formulas of the actual version of the coordinates can be found in Refs. 
. We use the same notations as Ref. [l(3|. For a grap hical view of the coordinates we 



refer the reader to Fig. 6 of Ref. 15| and Fig. 1 of Ref. 101] . In those studies the hyperangle 
changes in the interval [0; 2tt}. For systems with three identical particles (symmetry group 
Cs v ), which are only considered in this study, the possible irreducible representations are 
Ai, A 2 , and E. In order to represent A x or A 2 wave functions, it is enough to consider 
a smaller interval of 0, namely [7r/6;7r/2] (see Figs. 6 and 7 of Ref. [llj]). Wave func- 
tions of the E irreducible representation require a larger interval [vr/6; 7tt/6}. The systems 
considered here have Ai resonances only. However, the method works for other irreducible 
representations and can be adapted to systems with non-identical particles. The second 
hyperangle 9 changes in the interval [0;7r/2]. The shape of a three-particle configuration 
is determined by the two hyperangles, the overall size of the system is determined by the 
hyper-radius p G [0;oo). Therefore, the hyper-radius is the natural dissociation coordinate 
for the three-body system. 



B. The adiabatic hyperspherical approximation 

For a three-body system interacting through the potential V(p,9,<p), we have to solve 
the Schrodinger equation 

to obtain the vibrational three-body wave function $ n (p, 6, <p) and eigenenergy E^ b - In this 
study, V(p, 6, <fi) is given by a sum of two-body potentials. The present method does not 
use that property of the potential and works for potentials that are not separable into a 
simple sum of the three two-body terms. Usually, realistic three-body potentials cannot be 
represented in such a way. 
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$ n (p,e,<l>) = K b Mp,0,<f>) (2) 



One way of solving Eq. (j2j) is to use the hyperspherical adiabatic approach, treating the 
hyper-radius as an adiabatic coordinate. Then, Eq. (jSJ) is solved in a two-step procedure. 
In the first step, one 'freezes' the hyper-radius at pi and solves the Schrodinger equation in 
a two-dimensional space of hyperangles 

Hfoato e, 0) = u a { Pi ) Va { Pl - e, 0), (3) 

where <p a (pi', 9, 0) is the hyperangular wave function at hyper-radius p; for the a th eigenvalue 
U a (pi) (see, for example, Fig. [8]). The Hamiltonian H^f is obtained from the full Hamil- 
tonian of Eq. (j5J) by fixing the hyper-radius. The eigenenergies U a (pi) are called adiabatic 
potentials. Once the U a (pi) are obtained, Eq. (j2J) is written in the form of a system of 



coupled differential equations (see, for example, Ref. [201 ]) 



f(p) + U a (p)] Va,n(p) + WaAp)^a>M\ = K Va,n(p), (4) 

a' 

where T(p) is the kinetic energy operator associated with hyper-radial motion, W a , a '{p) are 
the coupling terms between the different channels a and a', and ip a ,n{p) is the a ih -channel 
component of the hyper-radial wave function. Assuming that the complete set of adiabatic 
channels a is included in Eq. (j4|), the solutions of Eqs. (j2j) and (HJ are equivalent. Sharp 
p-dependence of the coupling terms W a ^ a i (p) makes the numerical solution of Eq. (J3J) difficult. 

C. Mapped DVR method and slow variable discretization 

The slow variable discretization method suggested by Tolstikhin et al. \^ is used to 
account for the non- adiabatic coupling terms between the adiabatic states (p a {pi',9,4>) in a 
different way. We give only the formulas that are used in the current work. The detailed 
discussion of SVD can be found in Ref. 9|. The total vibrational eigenstate $(p, 9, 0) (index 
n is omitted here) is represented as an expansion in the basis formed by the adiabatic states 
tfiaiPi', 0) obtained from Eq. The expansion coefficients are ipa(Pi)'- 

$(p,0,<f>) = ^2il) a (pi)(p a (pi',6,<l>), (5) 

a 

where the sum is over all channels a. Expanding the hyper-radial wave function in a complete 
basis set ffj(p) 

Mp) = Yl c 3^Ap) i ( 6 ) 



we write the hyper-radial Schrodinger equation, Eq. (j3J), in the form of 

^ [(Tri\T{p)\iri>)Oia,i'a' + (■Ki\U a (p)\ir i >}5 aa > c Va > = E ^(71^)0^^^^ , (7) 

■V ,a' V ,a' 

where the overlap matrix elements Oi a ^ a i are 

O ia ,i<a> = (ipa{Pi\Q,4>)Wa'{Pi'\Q,4>)) ■ ( 8 ) 

These matrix elements account for the non-adiabatic couplings between the channels. They 
replace the p-dependent non-adiabatic terms W a>a i(p) in Eq. (j4j). The above Eq. (J7J) is 
written in a form of a generalized eigenvalue problem, which can be solved using standard 
numerical procedures. For an orthonormal DVR basis 7Tj, the equation is finally reduced to 

^2 [{n\T{p)\Tri')Oia,i'a' + U a {pi)5ii>5 aa > c^/ = E^20 ia>ia/ c ia t. (9) 



T(P) = --73- (10) 



In the above equations, the kinetic energy operator T(p) is 

fr^d, 2 
2p dp 

We use a DVR basis |7Tj) to represent the hyper-radial part of the wave functions. Since 
the weakly bound states and resonances have wave functions extending to large distances, 



we use a variable grid step for DVR basis functions similarly as we did in Refs. |10l . [16]. 
Namely, we introduce a new variable x such that p = p(x). The grid steps Ap and Ax along 
p and x, correspondingly, are linked to each other as Ap = J(x)Ax, where J(x) = dp/dx is 
the Jacobian. We may set Ax to be 1 and, therefore, Ap = J{x). In practice, first, we choose 
the dependence Ap(p). Once Ap(p) is chosen, the grid step in p automatically determines 
J(x). Then the integral J J(x)dx gives the dependence p(x). Changing the variable p to x 
in the kinetic energy operator T(p) and the wave function i[) to ip = y/Jipi Eq. becomes 

y^^Tu'Oia^g/Ci'a' + U a {p{Xi))c ia = E^2^ia,ia'Cia' , (H) 
i' ,a' a 1 

where Cj a = a/ J(xj)Q a . The matrix of the kinetic energy operator T(x) is determined once 
the DVR basis set 7Tj(x) is chosen. Previously, we used complex exponents, exp(ifcxj) as 



DVR basis functions. However, as it was pointed out in Refs. Ill, [l2|], the sine function 
basis, sin(ifcxj) has a certain advantage: zero boundary conditions at the ends of the grid are 
better represented by the sine functions than by the complex exponents. So, in this study 



we use the sine DVR basis. The actual form of the orthonormal DVR sine basis functions 
TTi(x) is well known. See, for example, Eqs. (5-7) of Ref. [11 1. To calculate the matrix 
elements Tm of the kinetic energy operator we use an approach similar to the one described 
in Refs. fill. \\\ . The matrix 7#/ can be represented as the following matrix product 

h 2 



T 



2// 



B A cos K A S i n RA sin K A cos B 



(12) 



where B, R, and K are diagonal matrices llj 



B; 



Sij/ J{xi) , Rij = 8ij/ J (xi) 



K -A — 



(13) 



N is the total number of the DVR basis functions. The elements of the matrices A S i n , A s ^, 



A C os, A^ s 



11 



12| 



.4. 



A s i n ) 



'(2J-1). 
sin | — — m | a,i , 



N 



sin 



(2j - 1; 

2N 



ITT OLi 



A, 



(2J-1). 
cos I in ] «j 

2 ( {2j-l) ( 

— cos tir a,- 

.V I 2.V 



(14) 



where 



1/V2, z = 0,AT, 
1 , otherwise , 
and i = 0, 1, . . .N;j = 1, . . .N . 



(15) 



To represent wave functions of loosely-bound states we choose a smaller step Api = J{xj) 
in the regions where the wave functions oscillate a lot, and a larger grid step for the regions 
where the wave functions change smoothly. One possibility is to connect the grid step Ap 
to the local kinetic energy E^ in (p) Q, 



Ap = /3- 



(16) 



^/2pE kin (p) 

where the parameter (3 < 1 provides additional flexibility. This choice for Ap gives good 
results in the regions where the WKB applicability condition is satisfied. In other situations, 
around classical turning points, for example, one needs to take a smaller step. 



D. Complex absorbing potential 



A complex absorbing potential is used to simulate an infinite hyper-radial grid (see Refs. 



171 ] . 181 ] and references therein). The CAP is placed at the end of the grid and, therefore, 
absorbs the outgoing dissociation flux of decaying resonant states. The length and the 
strength of CAP are chosen to minimize its effect: the outgoing flux should be completely 
absorbed and should not be reflected back to the internal region. The CAP is purely 
imaginary and added to the adiabatic potentials U a (p). This makes the Hamiltonian non- 
Hermitian. Solving the generalized complex eigenvalue problem of Eq. §7§ the eigenenergies 
of the Hamiltonian are obtained in the form 

E = E'-£, (17) 

where E' is the position of the resonance and V is the resonance linewidth (inverse of lifetime) 

Q. 

In literature, different types of CAPs are used (see Refs. [it]], [3] and references therein). 
In our calculations, we use quadratic and exponential forms of CAP. The widths and po- 
sitions of the resonances do not depend on the form of CAP if the length and strength 
are chosen appropriately (see Table 1171) . The optimal lengths and strengths of the CAPs 
are determined using the de Broglie wavelength of the dissociating products according to 



guidelines provided in Ref. 17J]. In this case, the dissociating products are a dimer plus a 
free atom. The kinetic energy of the dissociating products, from which we get the de Broglie 
wavelength, is the difference between the three-body bound state energy of the resonance, 
Etrimer, and the dissociation energy of the ground state adiabatic curve, Edi SS (see inset in 

Fig. m>. 

E. Reduced grid and variable grid-step in 2D space of hyperangles 

In this study, we are interested only in vibrational states of the A\ irreducible represen- 
tation. Such states can be represented by reduced interval [7r/6;7r/2] of the hyperangle (j> 
(examples of A\ functions are given in Figs. [H]and[H]). Using Neumann boundary conditions 

dy? a (0 = tt/6) _ <9y? a (0 = tt/2) _ Q ^ 
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we are able to represent states of the A\ irreducible representation. To represent states with 
A2 symmetry, we would have to impose Dirichlet boundary conditions. Such states come 
up in systems of fermions, for example (assuming that only the vibrational part of the total 
wave function determines the symmetry of the state). To represent E a or Ej, states, we 
would need to change our interval to [n/6,7n/2] and our boundary conditions according 
to E a (symmetric, Neumann) or Ej, (antisymmetric, Dirichlet) representation. 

We also employ a variable grid in hyperangles because the hyperangular wave functions 
are extremely localized when the hyper-radius is large. This occurs particularly in the 
asymptotic region of hyper-radius, where the hyperangular wave function is nonzero only 
in a small hyperangular region representing the dimer plus free atom configuration. For 
example, in the upper-left frame of Fig. [HI we see the hyperangular wave function very 
localized near 9 = 71/ 2 and = tt/6. The variable step size in the hyperangles, used in 
this study, is shown in Fig. [TJ The step size in 9, becomes smaller as 9 approaches 7r/2, 
while the step size in 0, approaches zero as approaches tt/6. In principle, it is possible to 
use different dependences like those shown in Fig. [T]for different hyper-radii. It would save 
some computation time. However, in this study, we use the same hyperangular grid for all 
values of pj. 




§ e 

FIG. 1: Variable grid step-sizes, A9 and A0, for 9 and 0, respectively. For a more accurate 
representation of hyperangular wave functions (see Figs. [U [9]), A0 — > as — ► 7r/6, whereas A 6 
-> as 9 -> vr/2. 
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III. EXAMPLE: THREE PARTICLES INTERACTING THROUGH A TWO- 
BODY POTENTIAL WITH A BARRIER 



A. Our results 

To test our method, first, we applied it to a simple model system that consists of three 
bosons with nucleon mass, m\ = = m 3 = 1837.5773 a.u. (939 MeV). The system was 



previously considered in Ref. [13j, where the Faddeev equations were employed in order 
to obtain three-body states. The three-body potential V{r\2, r23, r3i) in the Schrodinger 
equation, Eq. ([2]), is constructed as a sum of pairwise potentials: V^(r 12 , r 2 3, r 31 ) = V^r^) + 
^2(^23) + ^2( r 3i)) where the pairwise interaction V-^r) is given by: 

V 2 (r) = -55e- a2 '' 2 + 1.5e- a01 ^ 2 , (19) 

in MeV and the distance r is in fm. The two-body potential has a two-body bound state at 
£( 2 ) = -6.76 MeV Q. 

As the first step, we obtain the adiabatic energy curves. The resulting adiabatic potentials 
are shown in Fig. [3 The lowest potential at p — > 00 corresponds to a decay into a bound 
pair of two nucleons and a free nucleon. The dissociation energy of the ground state (Edi SS 
in Fig. [2]) corresponds to the dimer binding energy E™. All other channels dissociate into 
three free atoms. If the energy of a three-body state is above the lowest dissociation limit, 
the system may dissociate, with the excess of energy transferred into the kinetic energy of 
the dissociating products. 

With our method, we are able to calculate three-body resonances and bound states. We 
find that there is only one three-body bound state at En = -37.35 MeV. This result is in 

n 

agreement with the previous study by Fedorov et al. [13] giving En = -37.22 MeV. Figure [3] 
shows the multi-channel hyper-radial wave function ip{p) of this state. 

The pairwise interaction potential has a barrier. It means that if we have three particles 
situated close to each other, the dissociation to the dimer + free particle configuration is 
separated by a potential barrier too. Thus the lowest adiabatic state has a potenital barrier 
too. If the barrier is high enough, one or more resonances (pre-dissociated states) may exist 
with the energy above the lowest dissociation limit (dimer + free particle). The energy 
interval between the dissociation limit and the top of the barrier is roughly from -7 MeV to 
-5 MeV. 
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FIG. 2: (Color online) Hyperspherical adiabatic potentials for the system with a potential barrier. 
Five A\ adiabatic channels are displayed. There is only one bound state at -37.35 MeV. Resonant 
states are long-living states that tunnel through the potential barrier as shown in the inset. Note 
logarithmic scale in hyper-radius. The two arrows along hyper-radius point to the values of p where 
ground adiabatic curve has its minimum and to where CAP begins. 

Most of the eigenvalues of Eq. (TlT]) are neither bound states nor resonances: these are 
just states of the continuum spectrum. To distinguish the resonances from the continuum 
states, we may look into their wave functions or carry out calculations for a different grid 
length and CAP parameters of the grid: Positions and widths of resonances should be stable 
with respect to the change of the grid length and CAP parameters. In our calculations, we 
used CAP length and strength equal to 0.0028 a.u. and 80,000 a.u., respectively. We found 
only one resonance for this system with energy -5.31 MeV and halfwidth 0.12 MeV. The 
stability of this resonance with respect to varying CAP parameters is illustrated in Table 
HI We varied CAP length and strength parameters by +/ — 30% of optimal values. The 
wave function of the resonance is shown in Fig. HI The ground state component of the 
hyper-radial wave function has non-negligible components both inside the bound region as 
well as in the continuum region, as expected for a resonant state. Note that the ground state 
component of the hyper-radial wave function, ipi, is at least one order of magnitude greater 
than any other component, meaning that this state 'lives' mostly in the ground adiabatic 
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FIG. 3: (Color online) Ground state hyper-radial wave function, ip(p), broken down into compo- 
nents for five-channel calculation. Upper-left frame shows modulus squared of normalized wave 
function. Other three frames show channel by channel components, a = 1 is ground channel, a = 2 
is first excited state channel, and so on. Note logarithmic scale in hyper-radius for all hyper-radial 
wave functions. 

curve. Nevertheless, the presence of other components means that non-adiabatic couplings 
are not negligible. The significance of the non-adiabatic coupling is also demonstrated in 
Fig. [61 which shows the convergence of the results with respect to a different number of 
adiabatic channels included in Eq. (fTTl) : One channel is not enough to obtain an accurate 
result, two channels provide a much better accuracy. To compare with Fig. HI a continuum 
wave function for vibrational state with energy E = —6.4 MeV is shown in Fig. [5j An 
important difference with the resonant wave function is that the continuum wave function 
is situated mainly in the outer region of the hyper-radius, not reaching the potential well. 

For two-body systems, the complex energy of resonances does not depend on the shape of 
CAP, if CAP is smooth enough. So, we have also checked the convergence of our results with 
respect to different types of CAPs. Namely, we compared the quadratic and exponential 
CAPs. As expected, we found that the energy of the only resonance depends only weakly 
on the CAP parameters. Table [III compares the results obtained with the quadratic and 
exponential CAP. 



12 



"1 1 

- v 3 


1 1 11 





Ill 








(\ 1 1 1 












0.4 
0.2 


-0.2 
-0.4 



hyper-radius (a.u.) 



-5 ♦ -4 -3 

10 10 10 

hyper-radius (a.u.) 



FIG. 4: (Color online) Resonance hyper-radial wave function for four-channel calculation. The 
wave function has a considerable amplitude in the region of the potential well (around 10 -4 a.u.). 
Arrows along hyper-radius point to ground adiabatic curve minimum and to where CAP begins (see 
Fig. [2]). CAP is placed around p = 4 x 10~ 3 a.u., the distance at which the outgoing dissociation 
flux starts to be absorbed. This can be compared with the continuum wave function shown in Fig. 
[5l for which the potential-well amplitude is very small. Comparing the amplitudes of the different 
components, it is clear that the principal contribution is due to the lowest adiabatic state. 

B. Comparison with previous study 

Having developed the described method, we tested it on the present system because it 



was previously considered in Ref. [13J. We obtained the energy of the only bound state, 
which is in agreement with the mentioned study. However, the position (-5.319 MeV) and 
halfwidth (0.112 MeV) of the resonance do not agree with Ref. 13J, which gives position 
-5.96 MeV and halfwidth 0.40 MeV. This came as a surprise because the calculation for the 

n 

second three-body system considered here is in a good agreement with a previous work [8j . 
The reason of the disagreement is not clear. We are confident about the adiabatic states 
obtained in this study because we have obtained the same bound state as Fedorov et al. We 
are confident about the one-channel calculation giving a value of the complex energy, which 
is quite close to the converged value. The one-channel code does not use SVD at all and we 
tested it on diatomic molecules with a barrier. 
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CAP strength (a.u.) 


CAP length (a.u.) 


CAP initial position (a.u.) 


E'-i^ (MeV) 


80,000 


0.0028 


0.0060 


-5.307 - i0.116 


104,000 


0.0028 


0.0060 


-5.305 - i0.114 


56,000 


0.0028 


0.0060 


-5.310 - i0.117 


80,000 


0.0036 


0.0052 


-5.318 - i0.115 


80,000 


0.0020 


0.0068 


-5.344 - i0.107 



TABLE I: Stability of resonance with respect to CAP parameters for system of identical bosons 
interacting through a two-body potential with a barrier. The deBroglie wavelength of the dissoci- 
ating products is approximately 0.0006 a.u. The most deviant calculation, last row, is due to the 
short length of the CAP for that calculation. 



CAP type 


Energy (MeV) 


Halfwidth (MeV) 


quadratic 


-5.319 


0.1115 


exponential 


-5.318 


0.1183 



TABLE II: Stability of the calculated resonance with respect to type of CAP used. 

As an additional test, we carried out a calculation based on the WKB approximation. 
Strictly speaking, the WKB approximation can hardly be applied for this case, because 
the resonance is only the second state. WKB works only for excited states. But for the 
completeness of the discussion, we provide the details of the WKB calculation. In the WKB 
approximation, the linewidth is given by 



T 

The transmission probability Pt, and period T of oscillations are given by 

2 



In 



1 + exp 



T 



h 




2p,{U(p) - E trimer }dp 
dp, 



2p 



(20) 

(21) 
(22) 



Etrimer ~ U (p) 

where p is the three-body reduced mass, Etrimer is the energy of the resonance shown in Fig. 
([2]). For T, the limits of integration a and b are the classical turning points in the potential 
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FIG. 5: (Color online) Channel by channel components of continuum state wave function for 
four-channel calculation (E = —6.4 MeV.). Compared to the resonant wave function, Fig. [H 
the amplitude of this wave function inside the potential-well region is negligible. Arrows along 
hyper-radius point to ground adiabatic curve minimum and to where CAP begins (see Fig. [2]). 



well. The limits of integration b and c for Pt are the 'entrance' and 'exit' points through 
which the three-body dissociating flux is tunneling. In Fi g. E we compare our resonance 



linewidth to results previously obtained by Fedorov et al. 13j and WKB estimation. Our 



results are somewhat closer to the WKB approximation than those from Ref. 



131 ]. although 



as pointed out above WKB estimation does not provide a reliable result for this situation. 



One reason for the discrepancy with Ref. [13j could be that in our calculation we used a 
longer grid extending until p max ~ 0.01 a.u. than the grid in Ref. [TJ], where p max ~0.002 
a.u. We found that to obtain converged results within our approach, p max should be at least 
~ 0.005 a.u. or larger: CAP should have a length long enough to absorb the dissociation 
flux smoothly (see Figs. [3j HJ and [5]). In Ref. [131 ] . the situation could be different because 
no CAP was used. Instead, a complex scaling variable was employed. However, the complex 
scaling variable plays a role of an absorber. If the length of the complex scaling variable is 
too short, the outgoing dissociation flux is not absorbed completely and, thus, is reflected 
back from the end of the grid. This should give a wrong value for the resonance width. 

To provide an insight into the hyperangular part of the total wave function, in Figs. [8] and 
[9] we give hyperangular wave functions calculated for two different values of hyper-radius, 
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FIG. 6: (Color online) Convergence of the position and width of the resonance with respect to the 
number of adiabatic states included in the calculation. The results for the resonance are circled; 
other data points represent the continuum states. Positions and widths of the continuum states 
depend on the CAP parameters. Notice that positions of the resonance and continuum states are 
at negative energies because the origin of energy is at the three-body dissociation, which is above 
the dissociation to the dimer + free particle configuration. 

p = 0.001 and 0.01 a.u., correspondingly. For the small hyper-radius, the adiabatic states are 
delocalized. For the large hyper-radius, the lowest state is strongly localized: the amplitude 
is not zero only in the region where ~ 7r/6 and 9 ~ tt/2, corresponding to a dimer + 
free particle dissociation. Other adiabatic states at the large hyper-radius correspond to 
the three-body break-up: the corresponding wave functions are less delocalized, that means 
different possible rearrangements of the final dissociation products. 



a 
X 



0.01 




IV. EXAMPLE: EFIMOV RESONANCES IN A MODEL THREE-BODY SYSTEM 

The second system we consider is three interacting identical bosons that can form long- 
range quasi-bound states jsj], so-called Efimov resonances. The bosons have masses m = 1 
a.u. and interact through the two-body short-range potential 

V 2 (r) = -17.796 exp(-r 2 ). (23) 
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FIG. 7: (Color online) Comparison of our resonance linewidth calculations to linewidth calculation 
of Fedorov et al., and to the WKB approximation. We show our halfwidths for one-channel and 
five-channel calculations. Our four-channel calculation is closer to WKB approximation than that 
of Fedorov et al. [13]. However, WKB estimation cannot be viewed as reliable for this particular 
situation. 

(2) 

Two bosons with the interaction potential have a deeply bound state with energy E{ = 
—7.153 a.u. and a weakly bound state with energy Ef ] = -1.827 x 1(T 9 a.u. [8j. The 
s— wave scattering length is 23,394.87 a.u. The system of three bosons interacting 

through the potential with such a large scattering length and a relatively small radius of 
forces r ~ 1 is known to have a number of weakly bound three-body states jl|. Because 
for this particular potential there is the deeply bound two-body state, then weakly bound 
three-body states can decay into a dimer and a free particle. Therefore, the three-body 
states are indeed resonant states with finite lifetimes (and widths). 

Nielsen et al. [8] determined positions and widths of these states using an R— matrix type 
approach and the hyper-radial coordinates. Instead of solving the eigenvalue problem along 
the hyper-radius, in Ref. 8j the three-body scattering problem was considered. Namely, 
the energy- dependent scattering matrix S(E) was obtained for energies above the lowest 
dissociation limit (deeply bound dimer + a free particle). For the present situation when 
there is only one open channel, S(E) is a scalar function of E: S = e 2l&<yE \ The positions 
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FIG. 8: (Color online) Contour plots of hyperangular wave functions at p = 0.0001 a.u. for the 
system with potential given by Eq. (|19p . Values of contour lines are given by the legends. Upper 
left frame shows hyperangular axes: hyperangle 9 runs in the radial direction from to tt/2, while <j) 
runs in the polar angular direction from to 2ir. For this calculation we restict (f) to ir/6 < <p < 7r/2. 
This value of p corresponds to the bound state region. All wave functions are delocalized: there is 
no preferred three-body arrangement. 



and widths of the resonances were obtained considering the phase-shift 5(E). We use the 
results by Nielsen et al. [8j to test the present method. 

Since the Efimov states are very weakly bound, their wave functions extend to very large 
distances. For example, the minimum of the lowest adiabatic state is around 1 a.u., but the 
exponentially decaying tail in the closed channel of the fourth Efimov resonance reaches the 
distance of 10 5 a.u. (see Fig. 1 of Ref. 8j). At the same time, the component of the wave 
function corresponding to the open channel oscillates with the wavelength of A = 2.5 a.u. 
when hyper-radius is large. Therefore, in order to represent the complete wave function, one 
has to have a grid with a large number of points. Nielsen et al. js] used about 3 million grid 
points in the interval 0.1 < p < 10 6 a.u. 

Our method allows us to use a much smaller grid. The exponentially decaying tail can be 
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FIG. 9: (Color online) Contour plots of hyperangular wave functions at p = 0.001 a.u. for system 
with potential given by Eq. (|19h . Wave functions are represented in the same way as in Fig. El 
This value of p corresponds to the dissociation region. Note that the wave function of the lowest 
adiabatic state (left upper panel) is nonzero only in the small region of hyperangular space that 
represents the dimer plus free atom configuration for large hyper-radius. 



represented with a modest number of grid points (DVR points in our method). This can be 
made with a variable grid step, which grows logarithmically with hyper-radius. However, to 
be able to use the logarithmically-growing step, one has to absorb the outgoing dissociative 
flux at a small distance where the grid step is still small and the oscillating component of 
the wave function is still represented properly. Therefore, we place the absorbing potential 
at the lowest adiabatic potential only at a small distance. 

In a typical calculation, the variable grid in p has only 340 points and extends from 0.002 
to 20,000 a.u. The grid step is constant from p = 0.002 to 0.1 a.u. This region corresponds 
to the minimum of the two lowest adiabatic states (Fig. 1 of Ref. js]). From p = 1 a.u. 
the grid step changes logarithmically with p: Ap = 0.04p. We use a quadratic absorbing 
potential placed at the lowest adiabatic state. CAP starts at p =12 a.u; its length is L = 10 
a.u., and the strength A 2 = 11 a.u. The parameters L and A 2 are chosen according to Ref. 
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n 


E h units of 1(T 2 


E 2 , units of 10" 4 


£ 3 , units of 10~ 7 


E4, units of 10 9 


Ref. [s| 


-94.14 - i2.4705 


-28.38 - i3.001 


-58.42 - £5.79 


-19.92 - £1.771 


4 ad. states 


-96.0 - i2.38 


-28.8 - »2.91 


-56.7 -£5.59 


-22.6-il.67 


5 ad. states 


-96.0 - i2.39 


-29.1 - i2.98 


-57.2 -£5.72 


-22.5 - £1.76 


6 ad. states 


-96.0 - i2.39 


-29.2 - i3.00 


-56.3 -£5.75 


-22.6-il.77 



TABLE III: Comparison of complex energies (in atomic units) obtained in this work using different 
number of included adiabatic states with the results of Nielsen et al. 8j . The imaginary part of the 
energies is the halfwidth of the resonances. The overall agreement is good except for the position 
of the 4 th resonance, which is probably not represented accurately in this study: the grid is not 
long enough in the present calculation. In Ref. [8fl, six adiabatic states have been used. 



17l | in order to minimize reflection and transmission coefficients for the outgoing dissociation 
flux. In the calculation we used only 4 adiabatic states. The calculation of the four adiabatic 
states and the corresponding wave functions at 340 grid points of p is made in parallel on 96 
processors and took about an hour. The size of the matrix for the hyper-radial eigenvalue 
problem is (340*4) x (340*4). The construction and diagonalization of the matrix is made 
on a single processor and takes also about one hour. With the modest calculation effort, 
we obtained all four Efimov resonances with a reasonable accuracy. Table II III compares 
calculated energies and widths of the Efimov resonances with the results from Ref. jsj]. 
To give an idea about the wave function of an Efimov resonance, Fig. shows the four 
components of the third resonance. The wave function 'lives' mainly in the first excited 
adiabatic state: the probability to find the system in the first excited state is the largest. 
Regular oscillations are present mainly in the lowest adiabatic state, which corresponds to 
the open dissociation path. Beyond 20 a.u., the oscillations decay because of the complex 
absorbing potential. Smaller oscillations visible in other adiabatic states are due to the 
coupling to the lowest state. Oscillations along the states decay with the distance p because 
the coupling becomes smaller and the system can dissociate only through the lowest state. 
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hyper-radius (a.u.) hyper-radius (a.u.) 



FIG. 10: (Color online) First four components ip a (p) (a = 1 — 4) of the hyper-radial wave function 
of the third Efimov resonance. Real and complex parts of the components are shown in black and 
red lines. The main contribution to the wave function is from the second adiabatic tp a state having 
a minimum around 2 a.u. The only adiabatic state open for dissociation is the ground state, a = 1. 
However, the components with a = 2 — 4 have small oscillating tails (may not be visible), which are 
due to the coupling of the corresponding adiabatic states to ipi{p;9,4>). This is a generic property 
of the hyperspherical adiabatic approach: the couplings between adiabatic states decays slowly 
with p. Beyond p = 20 a.u., the oscillations on the lowest adiabatic state are damped due to the 
presence of the absorbing potential. 

V. SUMMARY AND CONCLUSIONS 

We developed a method to calculate three-body bound states and resonances. The 
method uses the adiabatic hyperspherical approach, slow variable discretization, sine DVR 
basis set, and a complex absorbing potential placed at a large hyper-radius. Using the 
hyper-radius as a dissociation coordinate allows us to account for two-body and three-body 
breakup uniformly. If needed, the three-body and two-body dissociation channels could 
be distinguished by inspecting the asymptotic behavior of resonance wave functions. To 
simplify the calculation, we employed a variable grid step in the two hyperangles and hyper- 
radius. In practice, the optimal variable change should be adapted to the given three-body 
potential. 
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We applied the developed method to two different three-body systems of identical parti- 
cles interacting through pairwise potential. However, the method works also for distinguish- 
able particles as well as for intrinsically three-body potentials. The first system consisted 
of particles interacting through the two-body potential with a barrier. Our calculation for 
the only bound state is in agreement with the previous study [13J , where the Faddeev equa- 
tions were used. However, the position and width of the only resonance does not agree 
with Ref. [ljj. The reason for this disagreement is not clear: we believe that our result is 
correct due to the larger size of the grid that we used. The second system considered here 
has four weakly bound resonances: Efimov resonances 8|. The widths and positions of the 
resonances obtained in the present study are in good agreement with results of Ref. . 

Among the advantages of the developed method is a relatively modest basis set along 
the hyper-radius. For example, in their .R-matrix method Nielsen et al. had to use several 
million grid points along the hyper-radius to obtain converged results. We needed only 
340 points. For a better precision, this number could be doubled. Another advantage is 
that the method does not rely on the evaluation of non-adiabatic coupling terms, which 
are difficult to represent accurately on a coarse hyper-radial grid. Most of the calculation 
can be done on a parallel computer, only the last step, which is relatively short, should 
be done on a single processor. This is an important advantage in view of the availability 
of parallel computers. The method works for intrinsically three-body potential and can be 
easily adapted for distinguishable particles. 

The method can be used in a number of problems dealing with long-living states of sys- 
tems with several degrees of freedom, if a dissociation coordinate could be separated from 
other coordinates. Such a dissociation coordinate should be used as an adiabatic coordi- 
nate. The other coordinates could be treated in the same way as the two hyperangles in 
the present approach. Depending on how good is the 'adiabaticity' of the dissociation coor- 
dinate, the number of adiabatic channels to include could be different: for the dissociation 
coordinate that is not adiabatic (slow), one may need to include many channels for conver- 
gence. Among possible applications of the method could be two-electron or electron-positron 
problems, pre-dissociation of small polyatomic molecules, decay of nuclei, processes in ultra- 
cold atomic gases. For example, applications to realistic atomic systems such as Efimov 
resonances in cesium p] will allow the exploration of the possibility of stabilizing the Efimov 
resonances into deeply bound levels by Raman processes. Such processes would allow the 
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formation of stable triatomic molecules. The present method provides accurate three-body 
wave functions, which are essential to compute the radiative transition probability. 
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